////////////////////////////////////////////////////////////////////////////////
//
// labFisPend.cc - acquire data for SNS Lab Fis pendulum experiment
// 
// 1st release: A. Lusiani <a.lusiani@sns.it> April 2003
//
// this program includes parts from poll.c, Part of Comedilib,
// Copyright (c) 1999,2000 David A. Schleef <ds@schleef.org>
// poll.c may be freely modified, distributed, and combined with
// other software, as long as proper attribution is given in the
// source code.
//

////////////////////////////////////////////////////////////////////////////////
//
// Modifications:
//
//

////////////////////////////////////////////////////////////////////////////////
//
// headers
//

extern "C" {
#include <stdio.h>
#include <stdlib.h>
#include <termios.h>
#include <comedilib.h>
#include <fcntl.h>
#include <unistd.h>
#include <errno.h>
#include <ctype.h>
#include <string.h>
#include <math.h>
#include <fcntl.h>
#include <sys/select.h>
}

#ifndef __CINT__
#include "TFile.h"
#include "TApplication.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TH2.h"
#include "TProfile.h"
#include "TAxis.h"
#include "TNtuple.h"
#include "TRandom.h"
#include "TStyle.h"
#endif

#include "labFisDaq.h"
#include "labFisPendEv.h"

////////////////////////////////////////////////////////////////////////////////
//
// function declarations
//

////////////////////////////////////////////////////////////////////////////////
//
// variable  declarations
//

const char* 		COMEDI_DEV = "/dev/comedi0";
//--- comedi sub-devices for Analog Input
int			AI_SUBDEV;

//--- buffer size for streaming data out of the driver
const unsigned int	COMEDI_BUFSIZ =	(4*1024*1024);

//--- Analog Input for light revealing photo-cell 
const unsigned int 	PHCELL_AI_CHAN=0;

//--- pendulum photo-cell analog input minimum/amaximum voltage
// const float 		PHCELL_VMIN = 0;
// const float 		PHCELL_VMAX = 5;
const float 		PHCELL_VMIN = 0.7;
const float 		PHCELL_VMAX = 2.3;

//--- Analog Input for temperatures
const unsigned int 	TEMP1_AI_CHAN=1;
const unsigned int 	TEMP2_AI_CHAN=2;

//--- Analog Input voltage range
const float 		AI_VMIN = 0;
const float 		AI_VMAX = 5;

//--- perform DAQ starting with AI #0, 3 channels
//--- #1 = photo-cell, #2 and #3 temperature probe
#define 		AISCAN_BEG	(0)
#define 		AISCAN_NUM	(3)
unsigned int 		aiScanList[AISCAN_NUM];

//--- select analog in voltage reference from GROUND
const unsigned int	AI_AREF = AREF_GROUND;

//--- rate at which the input is sampled, in Hz
unsigned int 		AISCAN_RATE = 50000;
#define 		AISCAN_SAMPLES (64*1024)

#define			HS_TIME_XAXIS	(200)

////////////////////////////////////////////////////////////////////////////////
//
// code
//

int main(int argc, char *argv[])
{
  int i;
  int rc;
  comedi_t *cifh;
  comedi_cmd cmd;
  lsampl_t maxDataAi = 0;
  double minVoltsAi, maxVoltsAi;
  unsigned int rangetAi;
  comedi_range* rangepAi = 0;
  //--- buffer for ni_pcimio Analog Input driver (16-bit unsigned samples)
  sampl_t buf[AISCAN_SAMPLES*AISCAN_NUM];
  lsampl_t phcellThreshUp;
  lsampl_t phcellThreshDown;
  
  //--- endpoints of acquisition are valid numbers
  rc = comedi_set_global_oor_behavior(COMEDI_OOR_NUMBER);

  cifh = comedi_open(COMEDI_DEV);
  if(!cifh){
    perror(COMEDI_DEV);
    exit(1);
  }

  //--- set comedi device streaming buffer size
  rc = comedi_set_buffer_size(cifh, AI_SUBDEV, COMEDI_BUFSIZ);
  if (rc<0) {
    fprintf(stderr, "comedi_set_buffer_size: %s\n",strerror(errno));
    exit(1);
  } else {
    rc = comedi_get_buffer_size(cifh, AI_SUBDEV);
    if (rc<0) {
      perror("comedi_get_buffer_size");
    } else {
      fprintf(stderr, "comedi buffer size for AI set to 0x%x bytes\n", rc);
    }
  }

  //--- find which subdevice offers Analog Input
  rc = comedi_find_subdevice_by_type(cifh, COMEDI_SUBD_AI, 0);
  if (rc<0){
    comedi_perror("comedi_find_device AI");
    exit(1);
  }
  AI_SUBDEV = rc;

  //--- find requested AI range, or most similar one
  rc = comedi_find_range(cifh, AI_SUBDEV,
			 PHCELL_AI_CHAN, UNIT_volt, AI_VMIN, AI_VMAX);
  if (rc<0){
    comedi_perror("comedi_find_range");
    exit(1);
  }
  rangetAi = rc;

  //--- get ptr to range struct
  rangepAi = comedi_get_range(cifh, AI_SUBDEV, PHCELL_AI_CHAN, rangetAi);
  //--- get maximum possible sample for AI
  maxDataAi = comedi_get_maxdata(cifh, AI_SUBDEV, PHCELL_AI_CHAN);
  minVoltsAi = comedi_to_phys(0, rangepAi, maxDataAi);
  maxVoltsAi = comedi_to_phys(maxDataAi, rangepAi, maxDataAi);

  //--- to detect up transition, ask threshold larger than 60% of range
  phcellThreshUp = comedi_from_phys(0.6*(PHCELL_VMIN+PHCELL_VMAX),
				    rangepAi, maxDataAi);

  //--- to detect down transition, ask threshold smaller than 40% of range
  phcellThreshDown = comedi_from_phys(0.4*(PHCELL_VMIN+PHCELL_VMAX),
				      rangepAi, maxDataAi);

  fprintf(stderr, "Comedi offers AI range n. %d, %f .. %f Volts\n",
	 rangetAi, minVoltsAi, maxVoltsAi);

  fprintf(stderr, "AI digital range: 0 .. 0x%x, threshold up 0x%x down 0x%x\n",
	  maxDataAi, phcellThreshUp, phcellThreshDown);

  //--- build channel list for DAQ scan of photo-cell
  for(i=0;i<AISCAN_NUM;i++){
    aiScanList[i] = CR_PACK(AISCAN_BEG+i, rangetAi, AI_AREF);
  }

  //
  // initialize Root
  //
  // code inspired from:
  // - test/hsimple.cxx of Root 3.05.3 package
  //

  TApplication theApp("App", &argc, argv);
  TCanvas *canvas = new TCanvas("canvas", "Pendulum Experiment", 600, 600);

  gStyle->SetTitleOffset(1.4,"Y");

  canvas->Divide(1,2);

  // Create a new ROOT binary machine independent file.
  // This file is now becoming the current directory.
  TFile hfile("labFisPend.root", "RECREATE", "Lab Fis: Pendolo");

  //--- allocate Root Tree
  TTree *tree = new TTree("nt", "Lab Fisica pendulum ntuple");

  //--- allocate pendulum event
  pendEv *event = new pendEv();

  //
  // define ntuple in Root tree
  // "ev" = name of Branch
  // "pendEv" = name of class used to store the event
  // 64000 = buffer size
  // 1 --> get a Root "leaf" for each member of class pendEv
  tree->Branch("ev", "pendEv", &event, 64000, 1);

  TH1F *hsPeriod = new TH1F("hsPeriod", "Oscillation period",
			   200, 1.95, 2.00);
  hsPeriod->SetXTitle("seconds");
  hsPeriod->SetYTitle("events");

  TProfile *hsShadowVsTime =
    new TProfile("hsShadowVsTime", "Shadow time vs. time",
		 HS_TIME_XAXIS, -HS_TIME_XAXIS/double(25), HS_TIME_XAXIS,
		 0, 2);
  hsShadowVsTime->SetXTitle("seconds");
  hsShadowVsTime->SetYTitle("seconds");
  // hsShadowVsTime->SetMaximum(0.4);
  // hsShadowVsTime->SetMinimum(0);
  hsShadowVsTime->SetMarkerSize(0.5);
  hsShadowVsTime->SetMarkerStyle(2);

  TProfile *hsPeriodVsTime =
    new TProfile("hsPeriodVsTime", "Period vs. time",
		 HS_TIME_XAXIS/10, -HS_TIME_XAXIS/double(25), HS_TIME_XAXIS,
		 0, 4);
  hsPeriodVsTime->SetXTitle("seconds");
  hsPeriodVsTime->SetYTitle("seconds");
  // hsPeriodVsTime->SetMaximum(1.990);
  // hsPeriodVsTime->SetMinimum(1.975);

  //--- prepare AI scan for microphone's input
  labFisScanPrepare(cifh, AI_SUBDEV, &cmd,
		    TRIG_NOW, AISCAN_RATE, TRIG_NONE, 0,
		    aiScanList, AISCAN_NUM);
  
  //--- ask for non-blocking I/O on comedi
  fcntl(comedi_fileno(cifh),F_SETFL,O_NONBLOCK);
  
  //--- allocate object to setup terminal input
  labFisTermIoHandler tioh(stdin);
  //--- setup terminal for immediate keypress input with no echo
  tioh.SetImmedNoecho();

  //--- start comedi Analog Input scan
  rc = comedi_command(cifh, &cmd);
  fprintf(stderr,"Start DAQ, rc= %d\n", rc);
  if (rc!=0) {
    fprintf(stderr,"errno=%d\n",errno);
    comedi_perror("comedi_command");
    exit(1);
  }

  canvas->cd(1);
  hsShadowVsTime->Draw("hist p");
  canvas->cd(2);
  hsPeriodVsTime->Draw("");
  canvas->cd();
  canvas->Update();
  
  int quit = 0;
  unsigned long counter = 0;
  unsigned long sampleCnt = 0;
  unsigned long counterLast = 0;
  unsigned long counterNew = 0;
  unsigned int counterOsc = 0;
  unsigned int periodOsc = 0;
  unsigned int shadowOsc = 0;
  unsigned int sample;
  unsigned int sampleLast = 0;
  unsigned int sampleNew = 0;
  unsigned int sampleNum = 0;
  unsigned int readMax = 0;
  unsigned int readMin = 0;

  unsigned int periodMin = static_cast<int>(1e7);
  unsigned int periodMax = 0;

  //--- DAQ loop
  while (1) {
    //--- if quit requested, end after completing cycle
    if (quit && (counterOsc %4 ) == 0) {
      //--- cancel endless DAQ
      rc = comedi_cancel(cifh, AI_SUBDEV);
      break;
    }

    //--- read data gathered by comedi DAQ
    rc = labFisScanRead(cifh, AI_SUBDEV, &buf[0],
			sizeof(buf),
			static_cast<long int>(1e6), 10);

    //--- update histograms
    canvas->Update();
    
    //--- return code == zero means end of data
    if (rc == 0) {
      fprintf(stderr, "EOF, ending scan after reading %6d transitions\n",
	      counterOsc);
    }

    //--- positive return code means number of available samples
    if (rc > 0) {
      unsigned int i;

      // fprintf(stderr, "read %8d samples, modulus %d\n",
      //	      rc/sizeof(sampl_t), rc/sizeof(sampl_t) % 3);

      //--- loop on scans, i.e. series of AISCAN_NUM samples
      //--- a partial scan (with n < AISCAN_NUM samples) might be left
      for(i=0;
	  i < (unsigned int)rc/sizeof(sampl_t);
	  i++, sampleCnt++, sampleLast = sampleNew, counterLast = counterNew
	  ){

	counter = sampleCnt/3;

	// fprintf(stderr, "0x%x %d\n", buf[i], counter);
	// fprintf(stderr, "%8d %x %x %x\n", i,
	//	buf[i], buf[i+1], buf[i+2]);

	if (sampleCnt % AISCAN_NUM ==1) {
	  //--- set event field with temp #1
	  event->t1 = comedi_to_phys(buf[i], rangepAi, maxDataAi);
	  continue;
	}

	if (sampleCnt % AISCAN_NUM ==2) {
	  //--- set event field with temp #2
	  event->t2 = comedi_to_phys(buf[i], rangepAi, maxDataAi);
	  continue;
	}

	//--- the following is only executed for photo-cell samples

	// fprintf(stderr, "%8ld %6d %x %x %x\n", sampleCnt, i,
	//	buf[i], buf[i+1], buf[i+2]);

	if (sampleLast) {
	  //--- detect up -> down transition
	  sample = (buf[i] >= phcellThreshDown);
	} else {
	  //--- detect down -> up transition
	  sample = (buf[i] > phcellThreshUp);
	}
	
	//--- first sample, initialize state machine
	if (sampleNum == 0) {
	  sampleNum++;
	  readMin = buf[i];
	  readMax = buf[i];
	  sampleNew = sample;
	  counterNew = counter;
	  continue;
	}
	
	if (buf[i] < readMin) {
	  readMin = buf[i];
	}
	if (buf[i] > readMax) {
	  readMax = buf[i];
	}
	
	//--- wait untill transition
	if (sample == sampleLast) continue;

	// fprintf(stderr, "%d %d %d %d\n", buf[i], buf[i+1], buf[i+2], counter);
	
	//--- note new sample and counter for each transition
	sampleNew = sample;
	counterNew = counter;
	sampleNum++;
	
	//--- discard 1st time difference (2nd sample)
	if (sampleNum <=2) continue;
	
	//--- begin acquisition with transition from dark to light
	if (sampleNum <=3 && sample > sampleLast) continue;
	
	if (sample < sampleLast) {
	  shadowOsc += counter-counterLast;
	}

	//--- fill event with light/ darkness periods
	switch (counterOsc % 4) {
	case 0:
	  event->l1 = counter-counterLast;
	  break;
	case 1:
	  event->o1 = counter-counterLast;
	  break;
	case 2:
	  event->l2 = counter-counterLast;
	  break;
	case 3:
	  event->o2 = counter-counterLast;
	  break;
	}
	//--- fill event with time counter
	event->time = counter;

	//--- the event is filled, store it on Root ntuple
	tree->Fill();
	
	periodOsc += counter-counterLast;
	//--- transition: save time diff. w.r.t. last sample counter
	printf("%6d ", static_cast<int>(counter-counterLast));
	//--- print 4 times per line
	if ( (counterOsc++ %4 ) == 3 ) {
#if 0
	  printf("min=%d max=%d mean=%d thresh=%d ",
		 readMin, readMax, (readMin+readMax)/2, phcellThresh);
	  printf("%6d\n", periodOsc);
	  printf("period = %6d shadow = %6d\n", periodOsc, shadowOsc);
#endif	  
	  if (periodMin > periodOsc) periodMin = periodOsc;
	  if (periodMax < periodOsc) periodMax = periodOsc;

	  {
	    double dTime = double(counter)/double(AISCAN_RATE);
	    double dPeriod = double(periodOsc)/double(AISCAN_RATE);
	    double dShadow = double(shadowOsc)/double(AISCAN_RATE);
	    double dPeriodMin = double(periodMin)/double(AISCAN_RATE);
	    double dPeriodMax = double(periodMax)/double(AISCAN_RATE);
	    
	    if (dTime > 0.90* hsPeriodVsTime->GetXaxis()->GetXmax()) {
	      hsPeriodVsTime->LabelsInflate();
	    }
	    hsPeriod->Fill(dPeriod);
	    hsPeriodVsTime->Fill(dTime, dPeriod);
	    
	    if (periodMin == periodMax) {
	      hsPeriodVsTime->SetMaximum(dPeriodMax+0.01);
	      hsPeriodVsTime->SetMinimum(dPeriodMin-0.01);
	    } else {
	      hsPeriodVsTime->SetMaximum(dPeriodMax
					 +0.1*(dPeriodMax-dPeriodMin));
	      hsPeriodVsTime->SetMinimum(dPeriodMin
					 -0.1*(dPeriodMax-dPeriodMin));
	    }
	    
	    if (dTime > 0.90* hsShadowVsTime->GetXaxis()->GetXmax()) {
	      hsShadowVsTime->LabelsInflate();
	    }
	    hsShadowVsTime->Fill(dTime, dShadow);
	    hsShadowVsTime->SetMaximum(1.25* dShadow);

	    canvas->cd(1);
	    hsShadowVsTime->Draw("hist p");
	    canvas->cd(2);
	    hsPeriodVsTime->Draw("");
	    canvas->cd();
	    canvas->Update();

#if 0
	    printf("temp1= %8.4f temp2= %8.4f\n",
		   event->t1, event->t2);
#endif
	    printf("time= %8.4f sh= %8.4f per=%8.4f\n",
		   dTime, dShadow, dPeriod);
	  }

	  shadowOsc = 0;
	  periodOsc = 0;
	  fflush(stdout);
	  //--- if quit requested, end after completing cycle
	  if (quit) {
	    break;
	  }
	}
      }

      continue;
    }
    
    if (-rc == LFDAQ_STDINRDY) {
      unsigned char ch = fgetc(stdin);
      if (ch == 25) {
	fprintf(stderr,"  --- <control-Y> received, will quit! --- \n");
	quit = 1;
      } else {
	fprintf(stderr,"  --- Press <control-Y> to quit --- \n");
      }
      continue;
    }

    if (-rc == LFDAQ_TIMEOUT) {
      fprintf(stderr,"Timeout waiting for data \n");
      break;
    }

    break;
  }

  fprintf(stderr, "Please select \"File->Quit Root\" on Root window to quit\n");

  //--- start root GUI: "Exit ROOT" will return to the next statement
  theApp.Run(kTRUE);

  //--- Save all objects in this file
  hfile.Write();

  //--- Close the file.
  hfile.Close();

  return 0;
}
